%% process_field
function process_field(fespace1, fespace2, U)
x1 = fespace1.node(:,1);
y1 = fespace1.node(:,2);
x2 = fespace2.node(:,1);
y2 = fespace2.node(:,2);
shift_u2 = fespace1.Nb;
shift_p = 2*fespace1.Nb;
U = full(U);
u1 = U(1:shift_u2);
u2 = U(shift_u2+1:shift_p);
p = U(shift_p+1:end);
% interpolation
xg = linspace(min(x1), max(x1), 64);
yg = linspace(min(y1), max(y1), 64);
[X, Y] = meshgrid(xg, yg);
U1 = griddata(x1, y1, u1, X, Y);
U2 = griddata(x1, y1, u2, X, Y);
P = griddata(x2, y2, p, X, Y);
% vector field
figure("WindowStyle", "docked");
velocityMagnitude = sqrt(U1.^2 + U2.^2);
pcolor(X, Y, velocityMagnitude);
shading interp;
colorbar;
hold on;
step = 2;
h = quiver(X(1:step:end, 1:step:end), ...
    Y(1:step:end, 1:step:end), ...
    U1(1:step:end, 1:step:end), ...
    U2(1:step:end, 1:step:end), 2);
h.Color = "w";
title("Velocity Field");
axis("equal", "tight");
hold off;
% streamlines
figure("WindowStyle", "docked");
streamslice(X, Y, U1, U2, 2);
title("Streamlines");
axis("equal", "tight");
% contours
figure("WindowStyle", "docked");
contourf(X, Y, P, 20);
colorbar;
title("Pressure Contours");
axis("equal", "tight");
end